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CONJUGATE GRADIENT TYPE METHODS 
FOR LINEAR SYSTEMS 

WITH COMPLEX SYMMETRIC COEFFICIENT MATRICES* 

ROLAND FREUNDf 

Abstract. Wc consider conjugate gradient type methods for the solution of large sparse linear 
system Ax = 6 with complex symmetric coefficient matrices A = A T . Such linear systems arise 
in important applications, such as the numerical solution of the complex Helmholtz equation. Fur- 
thermore, most complex non-Hermitian linear systems which occur in practice are actually complex 
symmetric. We investigate conjugate gradient type iterations which are based on a variant of the 
nonsymmetric Lanczos algorithm for complex symmetric matrices. In particular, we propose a new 
approach with iterates defined by a quasi-minimal residual property. The resulting algorithm presents 
several advantages over the standard biconjugate gradient method. We also include some remarks 
on the obvious approach to general complex linear systems by solving equivalent real linear systems 
for the real and imaginary parts of x. Finally, numerical experiments for linear systems arising from 
the complex Helmholtz equation are reported. 

Key words, complex symmetric matrices, nonsymmetric Lanczos algorithm, biconjugate gra- 
dients, minimal residual property 

AMS(MOS) subject classifications. 65F10, 65N20 

1. Introduction. Conjugate gradient type methods — used in combination 
with preconditioning — are among the most effective iterative procedures for solving 
large sparse nonsingular systems of linear equations 

(1.1) Ax = 6. 

The archetype of these schemes is the classical conjugate gradient algorithm (CG 
hereafter) of Hestenes and Stiefel [20] for Hermitian positive definite matrices A. 

While most linear systems which arise in practice have real coefficient matrices 
A and real right-hand sides fc, there are some important applications (see [14] and the 
references therein) which lead to complex linear systems. Partial differential equations 
which model dissipative processes usually involve complex coefficient functions and/or 
complex boundary conditions (see e.g. [23]), and discretizing them yields linear sys- 
tems with complex matrices A. A typical example for this category is the complex 
Helmholtz equation 

(1.2) -At i - cr x u + iVjti = /, 

where <T \ , <72 are real coefficient functions, which describes the propagation of damped 
time-harmonic waves as e.g. electromagnetic waves in conducting media. Further 
applications, which give rise to complex linear systems, include discretizations of 
the time-dependent Schrodinger equation using implicit difference schemes, inverse 


* This paper is based on a presentation at the Copper Mountain Conference on Iterative Meth- 
ods, April 1-5, 1990, Copper Mountain, Colorado. This work was supported in part by DARPA via 
Cooperative Agreement NCC 2-387 between NASA and the Universities Space Research Association 

(USRA). 

t RIACS, Mail Stop T045-1, NASA Ames Research Center, Moffett Field, CA 94035, USA, 
and Institut fiir Angewandte Mathematik und Statistik, Universitai Wuixburg, D-8700 Wurzburg, 
Federal Republic of Germany. 
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scattering problems, underwater acoustics, eddy current computations [2], numerical 
computations in quantum chromodynamics, and numerical conformal mapping. 

In all these examples, the resulting coefficient matrices A are non-Hermitian. 
However, they still exhibit special structures. Often, A differs from a Hermitian 
matrix only by a shift and a rotation: 

(1.3) A = e l *(T + id), T = T H Hermitian, <r 6 C, 9 6 R, : := V—T. 

In almost all other cases, which lead to complex systems, the coefficient matrix is 
symmetric: 

(1.4) A = A t is complex symmetric. 

Note that the two families (1.3) and (1.4) overlap. The matrix (1.3) is complex 
symmetric iff T is real. 

Surprisingly, when complex linear systems (1.1) are solved in practice, usually no 
attempt is made to exploit the special structures (1.3) or (1.4). Indeed, there are two 
popular approaches. The first one (see e.g. [1]) is to apply preconditioned CG to the 
Hermitian positive definite normal equations 

(1.5) A H Ax = A H b. 

Of course, complex numbers can always be avoided by rewriting (1.1) as a real linear 
system for the real and imaginary parts of x. The second popular approach is to solve 
this real and, in general, nonsymmetric linear system by one of the generalized CG 
methods such as GMRES [32]. It turns out that in both cases the resulting iterative 
schemes tend to converge slowly. As a consequence, complex linear systems have the 
bad reputation of being difficult to solve by CG type methods. On the other hand, for 
the class of shifted Hermitian matrices (1.3), there are efficient CG type algorithms [9, 
10, 14] for complex linear systems in their original form (1.1). We refer the reader to 
[14] for a detailed study and practical aspects of these schemes. In [14] it is also shown 
how the special structure (1.3) can be preserved by using polynomial preconditioning. 

In this paper, we are mainly concerned with CG type methods for linear systems 
(1.1) with coefficient matrices of the second class (1.4). In particular, we consider 
approaches based on a variant of the nonsymmetric Lanczos algorithm which was 
successfully used for computing eigenvalues of complex symmetric matrices (see [28] 
and [5, Chapter 6]). This Lanczos recursion generates basis vectors for the Krylov 
subspace induced by A which are orthogonal with respect to a certain indefinite inner 
product. The standard way to obtain from this basis iterates, which approximate the 
exact solution of (1.1), is to enforce a biconjugate gradient (BCG hereafter) condition. 
Here, we propose a new approach which generates iterates via a quasi-minimal residual 
property. On typical examples, the resulting algorithm displays better convergence 
properties than the BCG approach. In particular, it produces residuals whose norms 
are almost monotonically decreasing, in contrast to the erratic convergence behavior 
which is typical for BCG. Moreover, the new technique eliminates one of the two 
sources of possible breakdown in the BCG approach. 

The outline of the paper is as follows. In Section 2, we discuss the Lanczos re- 
cursion for complex symmetric matrices and state some of its theoretical properties. 
Section 3 deals with a variant of the BCG algorithm for the special matrices (1.4). 
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In Section 4, we propose the quasi-minimal residual approach for complex symmet- 
ric matrices. Section 5 contains some remarks on the problem of breakdown of t e 
complex symmetric Lanczos recursion. In Section 6, we are concerned with the issue 
“complex versus equivalent real linear systems”. In particular, some results are pre- 
sented which indicate that for Krylov subspace methods, such as CG type algorithms, 
it is always preferable to solve the original complex system rather than equivalent real 
ones. In Section 7, some typical results of numerical experiments for linear systems 
arising from finite difference approximations to the complex Helmholtz equation (1.2) 
are given. Finally, in Section 8, we make some concluding remarks. 

Throughout the paper, all vectors and matrices, unless stated otherwise, are 

assumed to be complex. As usual, M = [mjt], M T = [m*;]> and M H = M denote 
the complex conjugate, transpose, and Hermitian matrix, respectively , corresponding 
to the matrix M = [m ; *]. Moreover, we set ReAf = (M + M)/ 2 and ImM = 
(M — ~~Kf )/(2i) for its real and imaginary part, respectively. The notation 

K*(c, B) := span{c, Be, ... , B k ~ 1 c } 

is used for the Jfcth Krylov subspace of C n generated by c € C n and the n x n matrix 
B. Furthermore, c(B) denotes the spectrum of B. The vector norm ||x|| = vi x is 
always the Euclidean norm. The set of all complex polynomials of degree at most k 
is denoted by 

Hi := {p(A) = 7o+7i'* + --- + 7fcA ,! | 7o,7i> • • • .7* G C}. 

Moreover, the coefficient matrix A of (1.1) is always n x n and, unless stated otherwise, 
assumed to be complex symmetric. Generally, x* € C n , k = 0, 1, . • ., denote iterates 
for (1.1) with corresponding residual vectors r* := 6 — Ax*. If necessary, quantities 
of different algorithms will be distinguished by superscripts, e.g. xf 00 and x k 
Finally, an iterative scheme for solving (1.1) is called a Krylov subspace method if its 
iterates are of the form 

(1.6) x* e X 0 + Kk(ro,A) or, equivalently, x* = xo + P(A)r 0 , P G H*_i. 

Note that, in particular, CG type algorithms for (1.1) fall into this category. 

2. The Lanczos algorithm for complex symmetric matrices. In this sec- 
tion, we are concerned with a complex symmetric variant of the Lanczos algorithm 
and its theoretical properties. 

The Lanczos method [24] for general complex n x n matrices M is as follows (see 
e.g. [35, pp. 388-394]): 

Algorithm 2.1 (nonsymmetric Lanczos method). 

(1) Start: 

• Choose ro, So € C", ro, «o r ”* 

• Set v = ro, w = so, and vo = v>o = 0. 

(2) For Jb= 1,2,... do: 

• Compute T] = w T v; 

• // 1 } = 0: Stop; 

• Otherwise, choose fit , 71 G C with /?*7t = Vi 

• Set vt = v/ 7 * and wt = w/Pti 

• Compute at = wj Mvt; 

• Set r = Mvt — o*v* — fitvic-ii 

• Set w = M T wt - atm - 7kWk-i- 
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As Lanczos pointed out [26, p. 176], work and storage of Algorithm 2.1 can 
be halved if M is Hermitian resp. complex symmetric, by choosing starting vectors 
So = r 0 resp. so = ro. The resulting Hermitian Lanczos method has been studied 
extensively (see [18, Chapter 9] and the references therein). In contrast, the literature 
on the complex symmetric variant is scarce and restricted to the application of the 
algorithm to computing eigenvalues of complex symmetric matrices (see Moro and 
Freed [28] and Cullum and Willoughby [5, Chapter 6]). Here, we hope to convince the 
reader that the complex symmetric Lanczos algorithm is also very useful for solving 
linear systems. 

The basic method is as follows: 

Algorithm 2.2 (Lanczos method for A = A T ). 

(1) Start: 

• Choose ro € C n , ro 7^ 0; 

• Set v\ = ro and t>o = 0. 

(2) Fork = 1,2,... do: 

• Compute 0 k = (t5jv*) 1/2 ; 

• If P k = 0: Set mo = k — 1, and stop; 

• Otherwise , set v k = £>*//?*; 

• Compute Qk = vj Av^; 

• Set £*+1 = Av k - a k v k - 


In the next proposition, some elementary properties of Algorithm 2.2. are listed; 
proofs can be found in [5, Chapter 6]. We set 


( 2 . 1 ) 


V*:=[t>i v 2 ... «*] 
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Moreover, m* = m*(r 0 , A) := dim/C n (ro, A) denotes the grade of ro with respect to 
A (cl. [35, p. 37]). We remark that m* > 1 is the smallest integer such that is 
an A— invariant subspace of C n . Equivalently, if A is nonsingular and ro = b — Axo, 
m* > 1 is the smallest integer such that 


( 2 . 2 ) 


A 1 b E xo + K m ,(ro,A). 


Proposition 2.3. 

(a) In exact arithmetic, Algorithm 2.2 slops after a finite number of steps k = 
m 0 + 1 and 0 < m 0 < m*. Furthermore, t) mo+ i = 0 if mo = m* (“regular termina- 
tion”), and t5 mo +i ^ 0 if mo < m* (“breakdown”). 

(b) For Jfc = 1,2,..., mo: 


(2.3) 

(2.4) 

(2.5) 


T 

V k Vj 


ro, «/*#i 

l 1, ifk =j 


j = 1 , 2 ,... ,m 0) 


K k {r 0 ,A) = span{v 1( V2,... 

AV k = V k T k + [0 0 ... 0 wjt+i]. 
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Notice that, by (2.3-4), the Lanczos vectors form an orthonormal basis 

for /vt(r 0 , A) with respect to the (non-Hermitian) inner product 

(2.6) (x,y):=y T x, x,yeC". 

In particular, if Algorithm 2.2 terminates regularly, it generates a basis of the affine 
space xq + K m ,(ro , A) which, in view of (2.2), contains the exact solution of Ax = b. 

Next, we remark that (2.6) is the proper (cf. Craven [4]) inner product for 
complex symmetric matrices. Unfortunately, it has the defect that there exist vectors 
v £C n which are qua si-null [4], i.e. (v,v) = 0, but v ? 0. Consequently, it can not 
be excluded that Algorithm 2.2 actually breaks down. Indeed, in view of part (a) of 
Proposition 2.3, a breakdown occurs if one encounters a quasi-null vector v k . The 
phenomenon of breakdown will be discussed further in Section 5. 

We conclude this section with a result on the connection of the complex symmet- 
ric variant 2.2 with the general Lanczos Algorithm 2.1. Unlike Hermitian matrices, 
complex symmetric matrices do not have any special spectral properties. Indeed (see 
e.g. [21, Theorem 4.4.9]), any complex n x n matrix is similar to a complex sym- 
metric matrix. This result entails that the general nonsymmetric Lanczos Algorithm 
2.1 differs from the complex symmetric version 2.2 only in the additional starting 
vector so which can be chosen independently of t*o in 2.1. A strict statement of this 
correspondence is given in the following 

THEOREM 2.4. Lei M be a complex n x n matrix and ro € C n , ro ^ 0. 

(a) There exists a complex symmetric nxn matrix A which is similar to M: 

(2.7) M = XAX~ l where X is nonsingular. 

(b) Set r 0 = X _1 r 0 and s 0 = X~ T r a . Let v k ,w k ,a k ,0 k ,j k resp. v k ,a k ,p k be 
the quantities generated by Algorithm 2.1 (starting with r 0 ,so) resp. Algorithm 2.2 
(starting with r 0 ). Let m 0 denote the termination index for Algorithm 2.2. Then, for 
k = 1,2, .. . ,m 0 : 

(2.8) v k = (I| %-)x~ l v k = (n ^)x T w k , a k = at, (A)* = A 7*- 

Proof. Only part (b) remains to be proved. First, by means of (2.7), we rewrite 
Algorithm 2.1 in terms of A, X~ l v k , X T w k . By comparing the resulting iteration with 
Algorithm 2.2 and using induction on k, one readily verifies (2.8). 0 

After these preliminaries, we finally turn to linear systems (1.1) now. In the 
sequel, it is always assumed that A is nonsingular. 

3. The biconjugate gradient algorithm for complex symmetric matri- 
ces. In his celebrated papers [24, 25], Lanczos also proposed a scheme, closely related 
to Algorithm 2.1, for solving general non-Hermitian linear systems, namely the bicon- 
jugate gradient algorithm (BCG). We refer the reader to [11, 31, 22] for a detailed 
discussion of the BCG approach. 

Like Algorithm 2.1, BCG for general linear systems is started with two vectors: 
the residual r 0 = b - Ax 0 of the initial guess x 0 and a second vector s 0 t 0. We 
remark that s 0 is still unspecified. Due to the lack of a criterion for the choice of «o, 
one usually sets s 0 = r 0 in practice. For the case of complex symmetric matrices A, 
it is straightforward to show that, in analogy to the complex symmetric variant 2.2 
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of the general Lanczos Algorithm 2.1, the choice sq = ro results in a scheme which 
requires only half the work and storage of general BCG. The resulting procedure is as 
follows: 

Algorithm 3.1. (BCG for A = A r ) 

(1) Start: 

• Choose zq G C n ; 

• Set po = = b — Ax o and compute r^ro. 

(2) For k = 1,2,... do: 

• Compute Apk-i and pJ^iApk-i; 

• = 0 or rj_ l rj c ^i = 0: Set m\ = k — 1, and stop ; 

• Otherwise, set 6k = r *-i r *-i/pr_iApfc_i; 

• Compute x k = x*_i + £jfcP*-i and r* = r*_i - 6 k Ap k -i; 

• Comptife rJY* and set p k = r Jt / r Jl- 1 — 1 7 

• Compute p k = r* + pkPk-\- 

In the sequel, BCG always refers to the complex symmetric Algorithm 3.1. Next, 
we list some basic properties of BCG which follow readily from results (e.g. Jacobs 
[22]) for the general biconjugate gradient method. 

Proposition 3.2. 

(a) In exact arithmetic, Algorithm S.l stops after a finite number of steps k = 
mi + 1 and 0 < m\ < m*. Furthermore, x m , = A~ l b if mi = m* ( u regular termina- 
tion”), and x m t A *b if m i < m, (“breakdown’’). 

(b) Fork = 1,2 mi: 

(3.1) = 3 ~ 1.2, .- ,mi, 

(3.2) K k (r 0 , A) = span{r 0 , n, . . . , ft.,}. 

(c) Let k G {1,2, . . . T mi). Then , x k is uniquely determined by the Galerkin 
condition 

(3.3) (b-Ax k ) T y = 0 for all y€Kk(r 0 ,A), x k = x 0 + K k (r 0 ,A), 
with respect to the inner product (2.6). 

By comparing (3.1-2) with (2.3-4), we conclude that r*_i is parallel to the Lanc- 
zos vector v k generated by Algorithm 2.2. More precisely, one easily verifies that 

(3.4) r*_i = (— l)* -1 ii • • 0k-i0kVk, h = 1,2,... ,mi. 

Notice that there are two different causes for breakdown of Algorithm 3.1. The first 
one, namely the occurrence of a quasi-null residual vector r*_i, is, in view of (3.4), 
equivalent to the breakdown of the complex symmetric Lanczos Algorithm 2.2. In 
addition, Algorithm 3.1 breaks down if one encounters a search direction p k -\ ^ 0 
with Pk_ x Apk-\ = 0. This second cause of breakdown is more severe than the first 
one. As we will see in Section 4, it occurs if no Galerkin iterate (3.3) exists. 

Closely related to the biconjugate gradient method for general linear systems 
(1.1) is the conjugate gradients squared algorithm (CGS hereafter) which was recently 
proposed by Sonneveld [33]. 

Algorithm 3.3. ( CGS for general A) 

(1) Start: 

• Choose xo € C" and so G C n , «o ^ 0; 

• Set po = tio = r© = b — Azq and compute sj ro- 
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(2) For k= 1,2,... do: 

• Compute Apk-i and sq T A pt-i; 

• // $o T Apk-\ = 0 orsjn-x — 0; Stop; 

• Otherwise , set a* = sfi rk-i/$o Apu-i; 

• Compute qt = u*_i — cxkApk-ii 

• Compute Xk = x*-i + ajt(ujt-i + qk) and r* = r*_i — 

• Compute sjrjt and set (3k = 

• Compute ujt = rjt + (3t<lk; 

• Compute pk = u* + /?*(<?* + PkPk-x)- 


Notice that, like general BCG, CGS has a second unspecified starting vector s 0 . 
However, unlike BCG, even with the special choice s 0 = r 0j CGS can not exploit 
the complex symmetry of A. In particular, for A = A T , Algorithm 3.3 requires per 
iteration about twice as much work as the BCG Algorithm 3.1. 

Finally, as a special case of the general connection [33] between the CGS and 

BCG approaches, we have the following 

Proposition 3.4. Let A = A T , r 0 = r$°° = rtf™ , and , in Algorithm S.S, 

so = r 0 . Then, for fc = 0,1,... , mi, 

rf°° = Pk (A)r o and r«* = (p*(A)) 2 r 0 


/or some p* € 11* toi<A p*(0) = 1. 

4. A quasi-minimal residual algorithm for complex symmetric matri- 
ces. In this section, we propose a new approach for solving complex symmetric linear 
systems. The method is based on the complex symmetric Lanczos Algorithm 2.2. For 
simplicity, we assume throughout this section that, in exact arithmetic, Algorithm 2.2 
terminates regularly, i.e. 

(4.1) dt # 0 for 1.2,... ,m*, An«+i — 0- 

Moreover, let always be i € {1,2,... , m *} i n the following. 

4.1 Basic approach. Let x k be the Jkth iterate of any Krylov subspace method 
(1.6). Then, by (2.4) and with V k as defined in (2.1), we have 

(4 2) ** = *o + Viz* where z k € C*. 

Using (2.5) and r 0 = it follows from (4.2) that 

(4.3) r k =b- Ax k = r 0 - AV k z = ftt>i - V*+i f k z k — V*+i ifiiH - T k z k ). 

Here, t\ := [ 1 0 • • • 0] T denotes the first unit vector, 


(4.4) 



with ej := [0 • • • 0 1] , 


and, if Jk = m*, v m .+i := 0. Recall that 7* was defined in (2.1). 

Clearly, the aim is to choose z k in (4.2-3) such that r* « 0 as good as possible. 
In the BCG approach, this is attempted by enforcing the Galerkin condition (3.3). 
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Using (2.3-4) and (4.3), one easily verifies that (3.3) holds iff r* and v k +\ are parallel 
or, equivalently, z* is a solution of the linear system 

(4.5) TkZ = P\t\. 

Note that, by (4.1), (4.5) is inconsistent if T k is singular. Thus, we have the following 
Proposition 4.1. A BCG iterate zf 00 satisfying the Galerkin condition (3.3) 
exists if and only if T* is nonsingular . Moreover, if it exists, it is unique and given 

h 

(4.6) xf 00 = x 0 + V k Zk and r k = -/?*+i(z*)*v*+i 

where z k is the solution of (4.5) and (z*)* denotes its kth component . 

Proposition 4.1 demonstrates the defects of the BCG approach. Simple examples 
show that singular T* may indeed occur, and then, in view of Proposition 3.2, the BCG 
Algorithm 3.1 would break down in exact arithmetic. In floating-point arithmetic, such 
a breakdown is unlikely to happen. However, T* may still be close to singular and 
then the Galerkin condition (3.3) defines a poor approximation to the true solution 
of (1.1). This is the reason for the typical erratic convergence behavior with wildly 
oscillating residual norms. 

Obviously, the question arises whether there is a better strategy than (3.3) for 
choosing z* in (4.2-3). Ideally, one would like to have the minimal residual (MR) 
property 

(4.7) \\b-Axk\\= min ||6-Ax||= min ||V* + i (Aei - T*zfc)||. 

rexo+ZMro.A) z£C k 

However, by (2.3), in general (see Theorem 4.4 for an exception) the columns of V±+i 
are orthonormal only with respect to (2.6) rather than the Euclidean inner product 
(x,y) = y^x. Consequently, solving the least-squares problem on the right-hand side 
of (4.7) results in an algorithm for which work and storage per iteration step k grows 
linearly with k. Hence, if one insists on a “true” iterative scheme with constant work 
and storage per iteration, this excludes the MR method. 

Here, we propose the quasi-minima/ residua/ (QMR) approach as a substitute for 

(4.7) . Let 

= diag(wi, u> 2 , . . . , w*+i) with u>j > 0 for all j 
be a given positive diagonal weight matrix and rewrite (4.3) in the form 

(4.8) r t = [wifiiti -Slk+iTkZk)- 

Instead of ||r*|| as in (4.7), one may at least minimize the vector of components in the 
representation (4.8) of r*: 

(4.9) min Uwiftei - 

Hence, we define the iterates of the QMR method as follows: 

(4.10) x t = x < ^ tR = x 0 + V* 2 t where z* € C* is the solution of (4.9). 

Notice that, Jl * + \Tk is a (it + 1) x Jb matrix which, by (4.4) and (4.1), has full rank. 
Thus, the least squares problem (4.9) always has a unique solution z* . 

Clearly, the QMR approach still depends on the weights Uj . A natural choice is 

(4.11) Wj = |M, J = 1,2 k+ 1, 

so that all basis vectors Vj fuj in the representation (4.8) of r* have Euclidean length 1. 
Our numerical tests (cf. Section 7) also confirmed (4.11) as the best strategy. 
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4.2 Practical implementation. Next, we present an algorithm for the actual 
computation of the QMR iterates (4.10). The derivation is similar to that of Paige 
and Saunders’ SYMMLQ and MINRES algorithms [29] for real symmetric matrices. 

By (4.4), (2.1), and (4.1), flt+iT* is a tridiagonal (k + 1) x k matrix with full 
column rank. Hence, it admits a QR factorization of the type 

(4.12) = [ 0 ‘] 


where Qt+i is a unitary (k + 1) x (Jk + 1) matrix and R k a nonsingular matrix of the 
form 


(4.13) 


r Ci #3 


Rk ~ 


o C2 m 
0 C3 


0 


0 ■ 

0 

0* 


Lo 


• Vk 

0 <J 


The decomposition (4.12) can be generated by means of a series of k complex Givens 
rotations (e.g. [18, p. 47]) 


Q( c j> s j ) — 


Cj Sj 

. ~*j c i 


, Cj G R, Sj G C, Cj + |sj ! 2 — 1 » j — 1, • ■ ■ ,k. 


In particular, (4.12) is easily updated from the factorization QifliTt - 1 — R k -\ of the 
previous step by setting 


(4.14) 


n — f ^ 4-1 ® 

g * +1 -[o < 3(c*,st)J [ 0 lj 


and computing c k , s k and the new elements 6 k , ij k , Ct of Rk follows: 
0* = *lk — + 


(4.15) 


/ - 1/2 

£ k = C k -iU> k a k — S k -iC k -X*>k-l0k, |Ct| — (lCt| J + w *+llA+ll J > 

= r ictia/iai, if c* # 0 ct = < - t/Cti 8 k=Uk+l 0 k+ 1 /Ck. 

I iai. “ ‘ “ 


if Ct = 0, 

By (4.12) and since Q k +i is unitary, (4.9) is equivalent to 


min LAQ,+,.i - [**]*j- 


(4.16) 

From (4.10) and (4.16), it follows that 

(4.17) x k =x 0 + V k z k where z k := R k l t k , 


tk 

n+i 


| = tjk+l “ VlPlQk + l e l- 
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Notice that, in view of (4.14), tk differs from the previous vector tt- 1 only by its fcth 
component r t := (tk)k = ctf*. Next, we define vectors p, via 

(4.18) [pi P2 Vk]-=V k Rk l . 

Finally, using (4.17-18) and (4.13), one obtains the recursion 

Xk = x k -i + nPk, where pi = 

for the QMR iterates. In combination with Algorithm 2.2, the following implementa- 
tion results: 

Algorithm 4.2 (QMR method). 

(1) Start: 

• Choose xo € C " ; 

• Set vi — b — Ax o, vo = Po = P-i — 0; 

• Set 0i = (vf Cj)^ 2 , T\ = u>\0i , co = c_i = 1, and So — fi-i = 0. 

(2) For 4 = 1,2,... do: 

• If 0k = 0, slop: xi_i solves Ax = b. 

• Otherwise, compute Vk = i>k/0k ond a* = vjAvki 

• Sel v*+i = Av* - aivi - 0kVk-\, 0k+i = >' 

• Compute 6k, rji, Ci> c i> using formulas (4-15); 

• 5el pk = (t>* - rjkPk-i - 6kPk-i)/Ct; 

• Set Tk = ci f t , fi+i = -Skh; 

• Compute Xk — Xk-i + tupi. 


(Vl — »JlPl-.l - ^lPl-2) ! 


The assumption (4.1) guarantees that, in exact arithmetic, Algorithm 4.2 stops 
for k = m* + 1 and, by (2.2), n_i is indeed the solution of (1.1) then. However, 
in floating-point arithmetic, this finite termination property of the Lanczos recursion 
is no longer valid, and the stopping criterion stated in Algorithm 4.2 is not useful in 
practice. Instead, one should terminate the iteration as soon as ||r t || is sufficiently 
reduced. We remark that ri is not directly available in Algorithm 4.2. However, in 
view of (4.19), if one updates one additional auxiliary vector, namely 


hk = hk-i + 


c kh+i 

|SiS2---Sl| 2 U>l + l 


wt+i, 


h 0 •'= r 0 , 


then ||ri|| can be computed via 


INI = |sis 2 ---si| 2 -IIMI- 


Finally, notice that, for the weighting strategy (4.11), 


INI = 


\/s T s + t T t 

ifci 


s := Revt, 


Im i)i, 


can be obtained without extra cost during the computation of vj v* — s T s— t T t+2is T t . 
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4.3 Properties. In this subsection, we list some further properties of the QMR 
approach. 

Proposition 4.3. For k = 1,2,..., m*: 

(*) 

(4.19) rf MR = + (ctfk+i/wt+ijut+i; 

(b) The BCG iterate xf® 7 defined by (S.S) exists if, and only if, c* ^ 0. Moreover, 
if ct ^ 0, then 

(4.20) xf® 7 = + (h/ck)pk, 

(4.21) . Ikf 03 !! = • llwt+xllAwt+ilc*!). 

Proof, (a) By (4.17), (4.12), and (4.8), we have 


(4.22) 


r QM R _ f k+ iw k +i where w k +i ■= Vk+i^ZhQk+i 


0 

LI. 


With (4.14), it follows that successive vectors tDi+i and w k are connected by 


(4.23) Wk+i = — s k Wk + (ct/wi + i )«*+!■ 

Finally, by combining (4.23) and (4.22) and using fi+i = — sjtTt, one obtains (4.19). 
(b) First, we note that (4.12), (4.4), and (4.14) imply 


(4.24) 


QkftkTk = 


Ik-i 

0 



Rk. 


Thus, by Proposition 4.1, xf® 7 exists iff c* ^ 0. Now assume c t ^ 0. Using (4.5-6), 

(4.24) , and (4.17), we get 

(4.25) xf 03 = x 0 + Vizf® 7 where xf® 7 = R* 1 J - 

By comparing (4.25) with (4.17), (4.20) follows. For the proof of (4.21), notice that, 
by (4.25), (4.13) , and the formula for s k in (4.15), 

(4.26) fik+i(zk° G )k = 0k+ih/(CkC k ) = n**/(w*+ict). 


Furthermore, Algorithm 4.2 shows that 
(4.27) |?*| = |u>i0i«iS2 

Finally, by inserting (4.26-27) into the formula (4.6) for rf® 7 , we arrive at (4.21). □ 
In view of part (b) of Proposition 4.3, the QMB. method has the additional feature 
that it also yields a/I existing BCG iterates. This is in contrast to the BCG Algorithm 
3.1 which breaks down as soon as the first nonexisting BCG iterate is encountered. We 
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remark that, by (4.21), Hrf^H can be computed without extra cost from quantities 
which are generated in Algorithm 4.2 anyway. In particular, one may monitor ||rf ot? || 
during the course of the QMR algorithm, and, whenever the actual BCG iterate is 
desired, compute xf 05 via (4.20). 

There is an important special case for which the QMR method (with weighting 
strategy (4.11)) is even equivalent to the MR approach (4.7). Consider the subclass 
of (1.3) of complex symmetric matrices of the form 

(4.28) A = T+ i<r I, T -T r real symmetric, ff > 0. 

Assume that r 0 € R” (this can always be achieved by a proper choice of x 0 ). Then, it is 
easily verified that the Lanczos vectors v k generated by Algorithm 2.2 are all real and 
therefore, by (2.3), orthonormal with respect to the usual Euclidean inner product. 
In particular, by (4.11), w,- = 1, and the least squares problem (4.9) is equivalent to 
the one on the right-hand side of (4.7). Hence we have the following 

THEOREM 4.4. Lei A be of the form (4-28) an & r o G R n - Then, ihe iterates x k 
generated by Algorithm 4.2 (withw, = 1) satisfy ihe minimal residual property (4.7). 

5. On the breakdown of the complex symmetric Lanczos algorithm. Re- 
call that, throughout Section 4, possible breakdowns of the complex symmetric Lanc- 
zos recursion were explicitly excluded by assuming (4.1). In this section, we make 
some remarks about the general case and derive a theoretical result concerning so- 
called incurable breakdowns. 

First, let us return to the nonsymmetric Lanczos Algorithm 2.1. It stops as soon 
as w T v ~ 0 occurs. If this is caused by v = 0 or w = 0 , then one has found an invariant 
subspace. Unfortunately, Algorithm 2.1 may also break down, i.e. stop with w T v = 0 
and v, w ^ 0 (see e.g. [35, p. 389]). Although this happens very rarely in practice, the 
possibility of such breakdowns has brought the nonsymmetric Lanczos method into 
discredit and, certainly, kept many people from actually using the algorithm. However, 
especially due to the efforts of Taylor [34], Draux [6], Parlett, Taylor, and Liu [30], and, 
most recently, Gutknecht [19], the phenomenon of breakdown is now well understood. 
Moreover, there are look-ahead [34, 30, 19] variants of the Lanczos algorithm which 
allow to leap — except in the very special case of an incurable breakdown [34] — over 
those iterations in which the standard algorithm would break down. 

Here, we only sketch the basic idea of the look-ahead procedure for the special 
case of the complex symmetric Lanczos method. For a more detailed description of 
the look-ahead approach (for the general case) the reader is referred to [19]. Assume 
that breakdown occurs in Algorithm 2.2. In view of Proposition 2.3 this happens iff 
there is no complete set of m* Lanczos vectors v k € K k (ra,A), k = 1, — , m*, which 
are orthonormal (cf. (2.3)) with respect to the indefinite inner product (2.6). Clearly, 
there exists a maximal subset 

(5.1) {ki,k 7 ,...,kj}% {1,2,..., m*}, 1 < ij < k 7 < ■ • • < kj < m*. 

such that for each j =1,2,... ,J there exists a vector v kj € K k -(r 0 , A) satisfying the 
orthonormality relations 

(5.2) vj.v = 0 for all v 6 K kj -i(r 0 ,A) and vj.v kj = 1. 

By the definition of Krylov subspaces, K k (ro, A) = {P(A)r 0 | P G IU_i), and espe- 
cially 

(5.3) 


v kj = P kj -i{A)r 0 with P kj -\ € H^-j. 
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Therefore, we can rewrite (5.2) in terms of polynomials: 

(5.4) (Piy-i, P) = 0 for all P G IIij_2, P*,-i) # 0> 

with the indefinite inner product 

(5.5) (P,Q)--=r%P(A)Q{A)r 0 . 

A polynomial P kj . x € II t> _i which fulfills (5.4) is called a regular orthogonal (with 
respect to (5.5)) polynomial of degree kj - 1. It is well known [6, 19] that three 
successive regular orthogonal polynomials are connected via a three-term recurrence. 
By (5.3), it follows that there is a corresponding three-term recurrence relating the 
vectors v k -i, v t> , and v kj+l . The look-ahead Lanczos procedure is a modification of 
Algorithm 2.2 which — based on this three-term relation — generates the vectors v kj , 
j — 1,2, . . . , J ■ These vectors can then be completed to a basis of K kj by setting, e.g. 

v k = A k ~ k, v kj for k — kj + 1, kj + 2, . . . , ij+i — 1. j = 0, 1, . . • , J — 1- 

(cf. [17]). Here, for j = 0, we set k 0 := 1. We remark that the resulting look-ahead 
Lanczos algorithm produces block tridiagonal matrices T kj , j =1, . . . , J, of the type 

(2.1) with {kj - Jfcy _ i ) x {kj - kj. x) matrices a k] on the block diagonal. 

In exact arithmetic, the outlined algorithm terminates with the block tridiagonal 
T kj . Suppose that kj=m ic in (5.1). Then T kj represents the restriction of the matrix 
A to the A— invariant subspace I< m A r o, A). Obviously, in view of (2.2), the solution of 

(1.1) can then be computed from the quantities generated by the look-ahead Lanczos 
procedure. On the other hand, if kj < m* in (5.1), it is not possible to obtain the 
solution of (1.1) by means of the Lanczos process. For this reason, the case kj < m* 
is called incurable breakdown. 

Next, we derive a criterion for the occurrence of incurable breakdown in the 
complex symmetric Lanczos algorithm. In the following, it is assumed that A is 
diagonalizable. Then (e.g. [21, Theorem 4.4.13]), A has a complete set of orthonormal 
(with respect to (2.6)) eigenvectors. In particular, r 0 can be expanded into eigenvectors 
of A. More precisely, by collecting components corresponding to identical eigenvalues, 

we get 

••o = y\pw 

(5.6) i=i 

where pi ^ 0, Ati/ = Ajtcj, and, if l £ j> A; ^ A j, u t Uj — 0. 

Notice that, unless all eigenvalues of A are distinct, quasi-null vectors ui may occur 
in (5.7). In view of the following theorem, this is equivalent to incurable breakdown. 

Theorem 5.1. Let A = A T be a diagonalizable nxn matrix and r 0 € C". Then, 
in (5.1), kj = m, */, and only if, the eigenvectors in the expansion (5.6) of r 0 satisfy 

(5.7) uJui^bO for all 1= 1 m*- 

Proof. We need to show that (5.7) is equivalent to the existence of a regular 
orthogonal polynomial of degree m*- 1 with respect to the inner product (5.5). From 
the general theory of orthogonal polynomials, it is well known (e.g. [3, pp. 11 12J) 
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that regular orthogonal polynomials of degree k exist iff the corresponding moment 
matrix A/jt := ,* is nonsingular. For the case of (5.5), by (5.6), we have 


(5.8) fij := r%A 3 r 0 zz^pfXjuf ui, j = 0,1, 

1=1 

Note that moment matrices are in particular Hankel matrices. By applying Kro- 
necker’s Theorem on the rank of infinite Hankel matrices [16, pp. 204-207] to Afoo := 
(Mi+j)j,i=o, if follows that 

(5.9) rank A/qo — rank Af* = rank A/ m _i = m for all k > m — 1, 
where m is the number of poles of the rational function 

j=0 


Using (5.8) and Aj /x J ’ +1 = l/(z - A/), one obtains the following expansion of /: 

2 T 

(5.io) /(*) = ]P ~ r fora11 l 2 r l > I _. max l A, l- 


In particular, by (5.10), m < m* with equality holding iff (5.7) holds true. Hence, in 
view of (5.9), M m +-i is nonsingular iff (5.7) is fulfilled. This concludes the proof. □ 

As mentioned, (5.7) is guaranteed if A has only simple eigenvalues. Thus we have 
the following 

COROLLARY 5.2. If A — A T is an nxn matrix with n distinct eigenvalues , then 
incurable breakdowns can not occur in the complex symmetric Lanczos Algorithm 2.2. 

6. Complex versus equivalent real linear systems. In this section, we study 
connections between (1.1) and its equivalent real versions. Unless stated otherwise, A 
is now assumed to be a general complex n x n matrix. 

6.1. Equivalent real linear systems. By taking real and imaginary parts in 
(1.1), we can rewrite (1.1) as the real linear system 


( 6 . 1 ) 


A^ 


Rex 


Re6 

. [Re /I 

-ImyT 

Imx 

— 

Imi 

• ^ := [imA 

Rej4 


A second real version of (1.1) is 


( 6 . 2 ) 


Air * 


Rex 
— Imx 


Refrl A [Re A Im A 

Im6 * [Im A — Re A 


Obviously, (6.1) and (6.2) are the only essentially different possibilities of rewriting 
(1.1) as a real 2n x 2n system. Furthermore, note that A* is nonsymmetric iff A ^ A H 
is non-Hermitian, whereas A** is symmetric iff A = A T . Hence, for complex symmetric 
linear systems the approach (6.2) appears to be especially attractive since it permits 
the use of simple CG type methods such as SYMMLQ and MINRES [29] for real 
symmetric matrices. 



15 


In the following proposition, we collect some simple spectra) properties of A* and 

A**. 

Proposition 6.1. 

(a) Let J = X~ l AX be the Jordan normal form of A. Then A* has the Jordan 
normal form 


(6.3) 


J 0 
0 J 


= X~ 1 A ir X ir where A* := —r= 


v/2 l 


X - 

-iX 


-iX' 

X 


In particular, 

(6.4) <r(A+) = <r(A)Ucr(A). 


(b ) The matrices A** and —A** are similar. In particular, 

(6.5) —A, A,— A G <r(A**) for all A € <r(A**). 

Moreover, 

<r(A**) = {A € C | A 2 6 <r(AA)}. 

(c) Lei A - A t be complex symmetric. Then, there exists a singular value de- 
composition (the so-called Takagi SVD) of A of the form 

(6.6) A = ULU t , U unitary, E = diag(<ri, 02 , . . ■ , cr n ) > 0. 

Moreover, A** is a real symmetric matrix with spectral decomposition 


(6.7) 


Air* — 


y -z' 

E 

0 1 

y -z' 

Z Y 

0 

-eJ 

Z Y 


where Y = Re U, Z = Im U . 


Proof, (a) First, note that 


( 6 . 8 ) 



0 _ 

X 


where 


If in 

y/2 [“‘in 



is unitary. 


In particular, (6.8) 
that 


shows that with X also Af* is nonsingular. 

***-[• 3 ]' 


One readily verifies 


and, in view of (6.8), this implies (6.3). (6.4) is an obvious consequence of (6.3). 
(b) Since 



the real matrices A** and -A** are similar. Hence, (6.5) holds true. The relation 
between <r(A**) and <r(AA) is known (see [21, p. 214] for a proof). 

(c) (6.6) is the well-known Takagi singular value decomposition for symmetric 
matrices (e.g. [21, Corollary 4.4.4]). By rewriting (6.6) in terms of the real and 
imaginary parts of A and U , one obtains (6.7) (cf. [21, pp. 212—213]). D 

Roughly speaking, Krylov subspace methods are most effective for coefficient 
matrices A whose spectrum, except for possibly a few isolated eigenvalues, is contained 
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in a half-plane which excludes the origin of the complex plane. On the other hand, 
if this half-plane condition is not satisfied and if a large number of eigenvalues of A 
straddle the origin, usually the convergence of CG type algorithms is prohibitively 
slow. Typically, in these situations (see [7, 12, 13] for examples), iterations based on 
Krylov subspaces generated by A offer no advantage over solving the normal equations 
(1.5) by standard CG. See Theorem 6.4 for a theoretical result along these lines. 

For complex linear systems which arise in practice the half-plane condition is 
usually satisfied. Indeed, mostly 

(6.9) <r(A) C {A € C | ImA > 0}. 

However, by rewriting (1.1) as real linear systems (6.1) resp, (6.2), one deliberately 
creates coefficient matrices whose spectra are most unfavorable for Krylov subspace 
methods. The case (6.2) is especially bad since, in view of (6.5), <r(A**) is symmetric 
with respect to real and imaginary axis and hence the eigenvalues always embrace the 
origin. Similarly, by (6.4), the coefficient matrix A* of (6.1) in general has eigenvalues 
in the upper as well as in the lower half-plane. In particular, if (6.9) holds and, as 
in most applications, the Hermitian part (A 4* A**)f2 of A is indefinite, the spectrum 
of A* straddles the origin and the half-plane condition is not satisfied for A*. The 
following example illustrates this behavior. 

Example 6.2. Consider the class (4.28) of complex symmetric matrices A = T + 
i<xl where T = T T is real and c > 0. Obviously, 

ff(A) = {A = n + iff \ p G ff(T)} 

610 C S := [n m + iff, HM + «ff] • 

Here /i m and denote the smallest and largest eigenvalue of T, respectively. Note 
that the complex line segment 5 is parallel to the real axis and always contained in 
the upper half of the complex plane. In view of (6.4), (6.10) implies 

<t(A*) = {A = (i ± iff | fi € ff(T)} CSU5. 

We remark that Su2? is a tandem slit consisting of the two complex intervals 5 and S' 
which are parallel and symmetric to each other with respect to the real axis. Moreover, 
the eigenvalues of A* straddle the origin, if the Hermitian part T of A is indefinite. 
Finally, using (4.28) and part (b) of Proposition 6.1, we obtain 

ff(A**) = {A = ±vV + ff 2 | P G ff(T)} 

C [- y/fit U [ff.^/p^ + ff 2 j- 

Note that the class (4.28) is closely related to shifted skewsymmetric matrices. Indeed, 
if, instead of Ax — 6, we rewrite —iAx — —ib as a real system (6.1), one obtains 

(6.11) (-iA). = = ff/ 2 „ - N, N := Q ] (= -N 7 ^. 

Then, the eigenvalues are contained in a line segment which is parallel to the imaginary 
axis and symmetric with respect to the real axis: 

<r((-*A)*) = {A = ± i/i | \i e <r(T)} C [<r - *>, + ip), p = rnax{|/im|, Imat |}. 
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6.2. Correspondence of Krylov subspace methods. In analogy to (1.6) 
for complex linear systems (1.1), a Krylov subspace method for the solution of the 
equivalent real systems (6.1) resp. (6.2) generates iterates 


( 6 . 12 ) 


Rex* 

Imx* 


Rexo 

Imxo 


+ P(A*) 


Rer 0 
Imr 0 * 


Pen 


(r) 


resp. 

(6.13) 


Rex* 
— Imx* 


Re xo 
— Im xq 


+ P(A^) 


Rero 
Imro * 


Pen 


(r) 

Jb — 1 * 


Here and in the sequel, Ilj^j denotes the subset of II *- 1 of polynomials with real 
coefficients. Furthermore, the notation 


!<[%, B) := {P(B)c | P G (c K k (c, B)) 


will be used. 

At first glance, it might appear that Krylov subspace iterations (1.6) resp. (6.12- 
13) for the original complex systems resp. its equivalent real versions correspond to 
each other. However, as the following proposition shows this is not the case in general. 
Proposition 6.3. Let k e N. 

(a) Lei P E H*_i. Then, x* = x 0 + P(A)r Q is equivalent to 


(6.14) 


Rex*' 

Imx* 


Rex 0 

Imxo 




Rero 

Imro 




where P = Pi + iP 2 , Pj, P 2 € 11*2,. 

(b) Lei P € 11^2,. Then, (6. IS) is equivalent io 


(6.15) x* = Rex* + tlmn = xq + R(AA)tq + S(AA)Aro 


where R G an ^ & € n^_ 2 y 2 j aTt defined by P( A) = P(A 2 ) + AS(A 2 ). 

Proof. First, we note that, for j = 0, 1, . . 


(6.16) 



ReA' 
Im A> 


— Im A* 

ReX' 


and 



Re(yL4)> 
-Im (AAy 


Im (AAy 1 
Re (AAy J ’ 


as is easily verified by induction on j. 

(a) Let jj and 6, be the coefficients of the real polynomials P, and P 2 , respectively. 
Then, 


t-i 

ReP(A) = YXy Re A’ - 6j Im-A 2 ) 

i=o 

k - 1 

ImP(A) = YXfi A’ +6j Re A'). 

j=o 


(6.17) 


18 


By reformulating xt = xcH"-F(A)ro, by means of (6.17) and the first relation in (6.16), 
in terms of real and imaginary parts, one immediately obtains (6.14). 

(b) A routine calculation, using the second identity in (6.16), shows that (6.13) can 
be rewritten as 


Rex* 


Rezo 


— Imx* 


— Imzo_ 

T 


Re{R(AA)f 0 + S(AA)Ar 0 } 


Hence (6.13) and (6.15) are equivalent. □ 

In view of part (a) of Proposition 6.3, the corresponding real equivalent of com- 
plex Krylov schemes (1.6) are iterations of the type (6.14) and not the obvious real 
Krylov subspace methods (6.12). Clearly, the actual choice of the polynomials in (1.6) 
resp. (6.12-13) is aimed at obtaining iterates which are — in a certain sense — best 
possible approximations to the exact solution of the corresponding linear system. By 
using schemes of the type (6.12), from the first, one gives up k of the 2k real pa- 
rameters which are available for optimizing complex Krylov subspace methods (1.6.). 
Consequently, it is always preferable to solve the complex system (1.1) rather than the 
real version (6.1) by Krylov subspace methods. Furthermore, numerical tests reveal 
that the convergence behavior of the two approaches can be drastically different (see 
Section 7). 

6.3. A connection between MR and CGNR for complex symmetric 
matrices. Now assume that A is a complex symmetric nxn matrix. Then, in view 
of part (c) of Proposition 6. 1 , A** is a real symmetric indefinite matrix whose spectrum 
is given by 

(6.18) <r(A**) = {±*j | 3 = 1 **}• 

Here a, = <r,(A) > 0, j = 1, . . . , n, denote the singular values of A. 

Since there are simple extensions [29] of classical CG to real symmetric indefinite 
matrices, it is especially tempting to solve (6.2) by one of these methods. The iterates 
of these algorithms are defined either via a Galerkin condition or a minimal residual 
(MR) property. Here, we consider the MR approach. Applied to (6.2) it generates a 
sequence of iterates z*, k = 1,2 which are characterized by 

(6.19) 116** - A**z* || = min ||6** - A**z||, z* 6 *o + i£* r) (ro* , A**). 


Here, we have set 

(6.20) 6**:= [{£*], ** : = [-^mrj for * = := 6** - A**z 0 . 

Roughly speaking, CG type algorithms for real symmetric indefinite systems converge 
slowly if the coefficient matrix is strongly indefinite, in the sense that it has many pos- 
itive as well as many negative eigenvalues. Unfortunately, since, by (6.18), <r(A**) is 
even symmetric to the origin, A** exhibits this undesirable property. Indeed, numer- 
ical tests show that the convergence behavior of the MR method (6.19) is practically 
identical to that of the tabooed approach to (1.1) via solving the normal equations 
(1.5) by standard CG [20]. In the sequel, we refer to this latter method as CGNR. 
Notice that the iterates i* of CGNR are defined by the minimization property 

(6.21) 116 — Az/|| = min ||6-Ar||, x t € z 0 + Ki(A h r 0 , A 1 * A). 

v 9 " x£x 0 +Kt{A H r 0 f A H A) 
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Next, we prove that MR and CGNR are even equivalent, if the starting residual 
rj* satisfies a certain symmetry condition. Note that, corresponding to the spectral 
decomposition (6.7), rj* can be expanded into eigenvectors of A** as follows: 


( 6 . 22 ) 






' Ci “ 

Y 

-Z' 

c 

with c = 

; 

Z 

Y 



-C2n- 


Theorem 6.4. Let x^ R resp. xf GNR denote the iterates generated by (6.3.19- 
20) resp. (6.21) starting with the same initial guess xo € C" . Assume that c in the 
expansion (6.22) of r£* satisfies 


(6.23) ■ 

| c ; | — | c n+> |i i - 1,2,... ,n 

Then, 


(6.24) 

„OGNR _ -MR _ -MR / _ n i 

x x — x 7 1 —X 2 / + 1, « — v,i, 


Proof. First, note that, in view of (6.7) and (6.22), Cj and Cn+j are components 
corresponding to a pair of symmetric eigenvalues ±<r ; - of A**. However, for any real 
symmetric linear system A**z = 6** with “symmetric” eigenvalues and “symmetric” 
starting residual r$* in the sense of (6.18) and (6.23), respectively, the MR method 
generates iterates with zt G zo + (A**Cq* , Aj*) (see e.g. [13]). Consequently, 

the iterates defined by (6.19) satisfy 

(6.25) *2/ = *2f+i € *o + 

In particular, by (6.20), (6.25) shows that 

It remains to prove the first relation in (6.24). To this end, we remark that 

(6.26) ||6** - A**z|| = ||fc - A*|| for all z = [ , * € C B . 

Moreover, by using (6.20) and part (b) of Proposition 6.3 (applied to polynomials 
P( A) = AS(A 2 )), we deduce 

(6.27) z 0 + tf, (r) (A**rJ* , (A**) 2 ) = { [ | 1 € x ° + Kl r) (A H r 0 ,A H A)} 

(notice that A = A H in (6.15)!). In view of (6.25-27), (6.19) (for k = 2/) can be 
rewritten in the form 

(6.28) 116— Ax^H = min ||6-Ax||, a#* € x 0 +Kj r \A H r 0 , A H A). 

t€s,+X{ r) (A"r,,A"A) 

Finally, note that the iterates of CGNR always correspond to real polynomials, i.e. 
x oam € xo+kI t) (A h r 0 , A H A). Hence, by comparing (6.21) with (6.28), we conclude 
that if®” 1 = D 
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Clearly, the special symmetry condition (6.23) will not be satisfied in general. 
Nevertheless, all our numerical experiments showed (cf. Section 7) that (6.24) is still 
fulfilled approximately, i.e. 


(6.29) 


x 


CGNR „MR~MR 
; ~2 / ~2/+i> 


1 = 0 , 1 ,... 


7, Numerical Examples. We have performed numerical experiments with all 
algorithms considered in this paper in numerous cases. In this section, we present a 
few typical results of these experiments. 

Consider (1.2) on the unit square G := (0,1) x (0,1) with 0\ € R a constant 
and 02 a real coefficient function. First, assume that ti satisfies Dirichlet boundary 
conditions. Then, approximating (1.2) by finite differences on a uniform m x m grid 
with mesh size h := l/(m+ 1) yields a linear system (1*1) with A an n x n, n := m , 
matrix of the form 


(7.1) A-T + ihD, T := A o - 0 ih 2 I y D = diag(d x , d 2 , * * * >d n ). 

Here A 0 is the symmetric positive definite matrix arising from the usual five-point 
discretization of —A and the diagonal elements of D are just the values of 02 the 
grid points. 

Similarly, if we consider the real Helmholtz equation (1.2), i.e. 02 = 0, but now 
with a typical complex boundary condition such as 


— = iau on {(1, y) \ -1 < y < 1} 

and Dirichlet boundary conditions on the other three sides of the boundary of G, one 
again arrives at (7.1) where 


(7.2) 


f a, if j = /m, / sl,...,m, 
\ 0, otherwise. 


The test problems presented in this section are all linear systems Ax = b with 
complex symmetric coefficient matrices of the type (7.1). For Example 7.1, the mesh 
size h = 1/64 was chosen resulting in a 3969 x 3969 matrix A . In Examples 7.2-4, 
h = 1/32 and thus A is a 961 x 961 matrix. The right-hand side b was chosen to be a 
vector with random components in [—1,1]+ *[ — 1* 1]> with the exception of Example 
7.2 where b had constant components 1 + 1 . As starting vector xq = 0 was chosen. 

As stopping criterion, we used 


(7.3) 


_ 11* -Ax, 1| 

* ||6 — Ax 0 || 


< icr 6 . 


In Figures 7.1-4, the relative residual norm (7.3), R k , is plotted versus the iteration 
number i, at least for those methods for which work and storage per iteration is 
roughly the same. In the case of CGS resp. CGNR which both require about twice 
the work of the other algorithms and especially two matrix-vector products A * v resp. 
A • v, A • v per iteration, we have plotted R k versus 2 1. 

In a first series of experiments, QMR (with different weighting strategies) and 
BCG were compared. The natural choice (4.11) turned out to be the best strategy in all 
cases. In the following, QMR always refers to Algorithm 4.2 with weights (4.11). Then 
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QMR produces residual vectors whose norms are almost monotonically decreasing and 
generally smaller than those of the BCG residuals. However, convergence of QMR 
and BCG typically occurred after a comparable number of iterations. The following 
example is typical. 

Example 7.1. Here, (7.1) is a 3969 x 3969 matrix with <rj = 200, and the diagonal 
elements of D are given by (7.2) with a = 10. In Figure 7.1, the convergence behavior 
of BCG, QMR, and of the unweighted version (w> = 1) of the QMR Algorithm 4.2 is 
displayed. 


Figure 7.1 


Next, we compared the CGS Algorithm 3.3 with QMR and BCG. Typically, CGS 
needed slightly fewer iterations than QMR and BCG to reach (7.3). However, per 
iteration, QMR and BCG require only about half as much work and storage, CGS is 
not competitive for complex symmetric matrices. 

Example 7.2. In (7.1), we set n = 961, <n = 100 and dj,j = 1,-. . , n, are chosen 
as random numbers in [0, 10). Figure 7.2 shows the convergence behavior of QMR, 
BCG, and two runs of CGS with different starting vectors s 0 , namely s 0 = r o resp. 
So with random components in [—1,1] + »[— 1, 1). 


Figure 7.2 


Notice the extremely large residual norms in the early stage of the CGS iteration. 

In the following two examples, we compared CG type methods for Ax = 6 with 
real schemes for the equivalent real systems (6.1) resp. (6.2). For GMRES [32], work 
and storage per iteration step it grows linearly with k and in practice it is necessary 
to use restarts. In the sequel, GMRES(fco) refers to GMRES applied to (6.1) and 
restarted after every io iterations. Finally, MR(A**) denotes the minimal residual 
method (6.19) applied to the real symmetric system (6.2). 

Example 7.5. Here, in (7.1), n = 961, ffi = 100, and d s are given by (7.2) with 
a = 100. In Figure 7.3, for QMR, MR(A**), GMRES(5) resp. CGNR, the relative 
residual norm (7.3) is plotted versus iteration number k resp. 2 k. 


Figure 7.3 


Notice that, although the symmetry condition (6.23) is not fulfilled, the curves for 
CGNR and MR(A**) are almost identical. This confirms (6.29). Finally, we tried 
GMRES(Jto) also with other restart parameters to- For this example, the method did 

never converge. . 

Example 7 . 4 . Let A be the 961 x 961 matrix (7.1) with <r\ = 1000, D - ajl, 
<t 7 = 100 and set a := Note that A is a shifted Hermitian matrix of the form 

(4.28) (cf. Example 6.2). In particular, A belongs to the class of matrices (1.3) for 
which efficient true minimal residual algorithms for solving Ax = b exist. Here we used 
the particular implementation, MR(A), derived in [14, Algorithm 2]. Recall that, by 
rewriting — iAx — —ib as a real system (6.1), one obtains a shifted skewsymmetric 
matrix (6.11), (-L4)*. Again, for such matrices an efficient true minimal residual 
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algorithm, denoted by MR((-»j 4)*), exists [8, 12]. Figure 7.4 shows the convergence 
behavior of MR(^4), MR(A**), MR((— »j 4)*),CGNR, and GMRES(20). 


Figure 7.4 


Notice that MR((— and CGNR are nearly identical. This is typical for the case 
that a is small compared to the spectral radius of T. Furthermore, if a = 0, i.e. 
(~iAU in (6.11) is skewsymmetric, CGNR and MR((— :A)*) are even equivalent [12]. 

8. Concluding remarks. Complex linear systems Ax = 6 which arise in prac- 
tice often have complex symmetric coefficient matrices A . In this paper, we have 
explored the use of a variant of the nonsymmetric Lanczos process for complex sym- 
metrise matrices for the solution of such linear systems. In particular, we have proposed 
a new method of defining approximate solutions of Ax = b via a quasi-minimal resid- 
ual (QMR) property. In contrast to the biconjugate gradient (BCG) approach, the 
QMR iterates are well-defined as long as the basic Lanczos recursion does not break 
down. Moreover, unlike the wildly oscillating BCG residuals, the QMR residuals con- 
verge almost monotonically. Also, existing BCG iterates can be easily computed from 
the quantities generated during the QMR iteration. Finally, possible breakdowns — 
except incurable ones — of the complex symmetric Lanczos recursion can be overcome 
by using a look-ahead version of the Lanczos process. Incurable breakdowns only oc- 
cur in very special situations. For example, they can not occur if all eigenvalues of A 
are distinct. 

It is very tempting (and often done in practice!) to avoid complex linear system 
by solving equivalent real systems instead. We have presented some theoretical and 
numerical results which show that this — at least for Krylov subspace methods — is 
a fatal approach. Typically, the resulting real systems are unequally harder to solve 
by conjugate gradient type algorithms than the original complex ones. 

In this paper, we have not addressed the question of how to choose precondition- 
ers M for complex symmetric linear systems. This will be the subject of a forthcoming 
report. Here, we only remark that complex symmetry is preserved under precondi- 
tioning as long as M is complex symmetric. In particular, all algorithms for A = A T 
which we have considered can be used in conjunction with a complex symmetric pre- 
conditioner M . Note that the standard techniques, such as incomplete factorization 
[27], applied to A = A T generate complex symmetric preconditioners M. 

Finally, we would like to mention that the quasi-minimal residual approach can 
also be used to stabilize the general nonsymmetric biconjugate gradient algorithm [15]. 
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